Geometric frustration in small colloidal clusters 
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Abstract. We study the structure of clusters in a model colloidal system with 
competing interactions using Brownian dynamics simulations. A short-ranged 
attraction drives clustering, while a weak, long-ranged repulsion is used to model 
electrostatic charging in experimental systems. The former is treated with a short- 
ranged Morse attractive interaction, the latter with a repulsive Yukawa interaction. 
We consider the yield of clusters of specific structure as a function of the strength 
of the interactions, for clusters with m = 3,4,5,6,7,10 and 13 colloids. At sufficient 
strengths of the attractive interaction (around 10 fcsT), the average bond lifetime 
approaches the simulation timescale and the system becomes nonergodic. For small 
clusters m ^ 5 where geometric frustration is not relevant, despite nonergodicity, for 
sufficient strengths of the attractive interaction the yield of clusters which maximise the 
number of bonds approaches 100%. However for m — 7 and higher, in the nonergodic 
regime we find a lower yield of these structures where we argue geometric frustration 
plays a significant role, m = 6 is a special case, where two structures, of octahedral and 
C2v symmetry compete, with the latter being favoured by entropic contributions in 
the ergodic regime and by kinetic trapping in the nonergodic regime. We believe that 
our results should be valid as far as the one-component description of the interaction 
potential is valid. A system with competing electrostatic repulsions and van der Waals 
attractions may be such an example. However, in some cases, the one-component 
description of the interaction potential may not be appropriate. 



PACS numbers: 82.70.Dd; 82.70.Gg; 64.75. H-g; 64.60.My 
1. Introduction 

Clusters are a distinct state of matter which exhibit different structural ordering relative 
to bulk materials [1]. Of particular relevance to, for example, many biological systems 
such as viruses, is their tendency to exhibit five-fold symmetry such as icosahedra and 
decahedra [2] . Recently there has been a surge of interest in clusters formed in colloidal 



Geometric frustration in small colloidal clusters 



2 



systems Plll5l[6l[7llHlinilinilIIl[ll[15l[T8], which is expected to lead to the development 
of 'colloidal molecules' [U [191 ttOl |20l [15] . These may in turn provide novel functionalised 
materials d UHl [IDl [IH [2Q1 US] . 

Part of the attraction of studying colloidal dispersions is that, although in principle 
they are rather complex multicomponent systems, the spatial and dynamic asymmetry 
between the colloidal particles (10 nm-1 /xm) and smaller molecular and ionic species has 
led to schemes where the smaller components are formally integrated out |2lj. This leads 
to a one-component picture, where only the effective colloid-colloid interactions need be 
considered. This is usually a good approximation, however, there are a few exceptions. 
To describe the effective interactions between charged colloidal particles, the screening of 
the counter-ions is treated on a linear response level resulting in the well-known screened 
Coulomb pair potential proposed by Derjaguin, Landau, Verwey and Overbeek (DLVO). 
In general, however, due to non-linear counter-ion screening, the effective interactions 
involve many-body contributions [22], which violates the additivity of the potential. 
Furthermore, in strongly driven systems hydrodynamic interactions (velocity fields) 
modify the electrostatic interaction significantly [23| [2l] . Here we assume the one- 
component description. 

The colloid behaviour in the original complex system may then be faithfully 
reproduced by appealing to liquid state theory [2S] and computer simulation [2S]- 
Since the shape of the particles is typically spherical, and the effective colloid-colloid 
interactions may be tuned, it is often possible to use models of simple liquids to 
accurately describe colloidal dispersions. For example this approach has made it possible 
to reproduce Bernal spirals seen experimentally [HI [27] in a system with competing short- 
ranged attractions and long-ranged repulsions, leading to the idea that colloidal gels can 
be stablised by repulsive interactions [281 EHl . In addition to their own fascinating 
behaviour such as novel diffusion [T3] colloidal clusters are also predicted from theory 
and simulation to exhibit hierarchical ordering such as lamellae [31], cluster crystals [32] 
and may also undergo dynamical arrest to form cluster glasses [281 [SSI [33] . 

Since colloids may be directly imaged at the single-particle level, one may consider 
local intra-cluster structure, along with cluster-cluster correlations, a level of detail 
seldom accessible to atomic and molecular systems, except in the low-temperature 
regime [3l]. Meanwhile, the behaviour of colloidal clusters, for example the global 
energy minimum structure, should exhibit similarities to that of clusters of Noble gas 
atoms as both colloids and Noble gases can be well described by spherically symmetric 
attractive interactions [351 [36] • Some of us recently compared direct imaging of colloidal 
clusters with expectations from theory [37]. In the experiments, a significant deviation 
from expectations was found, in particular a maximum yield of only around 10% in 
structures expected to minimise the potential energy for small 4 ^ m ^ 6 clusters. This 
is surprising, as at these small sizes, there is little geometrical frustration that might 
inhibit access to the ground state, as would occur for larger clusters. 

It is the purpose of this work to establish those cluster structures we expect in the 
case of model colloidal systems by applying tried and tested model interactions. We 
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shall consider weak electrostatic repulsions and a short-ranged attraction for similar 
parameters to the experimental system. Typically, the former stems from electrostatic 
charge, the latter from the addition of non-adsorbing polymer or van der Waals 
attractions, and drives clustering. Rather than seeking minimum energy states [351 ES] , 
instead we shall try to mimic experimental approaches, as a way to predict what sort 
of yields of desired clusters we find as a function of interaction strength. Experiments 
on colloidal systems typically start from a randomised initial state. Unlike molecular 
systems, the effective temperature is typically constant: temperature itself is not 
usually varied and the effective colloid-colloid interactions in a given sample are taken 
to be fixed, as they depend upon sample composition. In other words the system 
undergoes an 'instantaneous quench' from an initial, randomised state. We follow this 
protocol, although we note recent theoretical and simulation work highlighting the role 
of microscopic reversibility in optimising yields in self- assembling systems |38l [391 SO] 
which is beginning to be exploited in nanoscience [H]. For the larger cluster sizes 
considered, we also investigate the role of electrostatic repulsions in cluster elongation, 
as suggested theoretically |12] and noted experimentally |5l |6l [7] . 

Our approach is as follows. We shall consider clusters of m = 3,4,5,6,7,10,13 
particles and study the structures formed, with particular reference to the minimum 
energy states ES], as a function of the attractive interaction strength (which 

mimics the addition of polymer in experimental systems). Various strengths of the 
repulsive interaction (which models the electrostatic charge in experimental systems) 
are also considered. For high strengths of the attractive interaction (low effective 
temperature), one expects the majority of clusters to reside in their minimum energy 
state. However, average bond lifetimes can exceed the simulation run time (which is 
comparable to experimental timescales) for sufficiently strong interactions. The system 
is thus nonergodic on these timescales, and kinetic trapping may become important. For 
larger clusters, long-ranged repulsions may lead to elongation [HI El [71 [271 133 132] • We 
investigate this effect for m = 10, 13. Our main result is that for small clusters (m ^ 5) 
the lack of geometric frustration enables the minimum energy state to be accessed, while 
for larger clusters kinetic trapping is important. 

This paper is organised as follows: section [2] introduces the model interactions, our 
simulation methodology is presented in section [3l followed by results (section [1]) and a 
discussion in section |5] in which we place our findings in the context of recent work. We 
conclude our findings in section [6l 

2. Model 

The seminal theory of colloid-polymer mixtures is that of Asakura and Oosawa [13] . 
This AO model ascribes an effective pair interaction between two colloidal hard spheres 
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in a solution of ideal polymers which is plotted in figure [H and reads, 
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where /3 = l/ZcgT, r is the centre to centre separation of the two colloids and the 
polymer fugacity zpr is equal to the number density ppp of ideal polymers in a 
reservoir at the same chemical potential as the colloid-polymer mixture. Thus within the 
AO model the effective temperature is inversely proportional to the polymer reservoir 
concentration. The polymer-colloid size ratio q = 2RG/cr where Rg is taken as the 
polymer radius of gyration, and a is the colloid diameter. For small polymer-colloid size 
ratios, the AO model has been found to give good agreement with direct experimental 
measurement [H]. 

The discontinuous nature of the AO interaction at contact complicates its use in 
Brownian dynamics simulations. We therefore use the continuous Morse potential, 
which, for short interaction ranges, is very similar to the Asakura-Oosawa potential 
(figure [1]) |15]. The Morse potential reads 



where po is a range parameter and /Sem is the potential well depth. The global energy 
minimum structures for clusters interacting via the Morse interaction are known [35] , 
and for small clusters m < 8, the structure of the global energy minimum is not sensitive 
to the range of the interaction. With the exception of m = 6 (see below), small ground 
state clusters should be the same for an AO colloid-polymer system as those tabulated 
for the Morse interaction. The experimental system [37j used a polymer-colloid size 
ratio of q = 0.22 which maps to a Morse range parameter po = 33.06 for a well depth 
(3eM = 5.0 according to the extended law of corresponding states [55]. 

Repulsions in colloidal systems typically stem from the electrostatic charge on the 
colloidal particles. Under many conditions, where the charging is quite weak, as is the 
case here, this leads to a screened Coulomb, or Yukawa form with a hard core that 
accounts for the physical size of the colloidal particles: 



where k is the inverse Debye screening length. The contact potential is given by 
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Figure 1. (color online) (a) The Morse potential [eq. ^] matched to the Asakura- 
Oosawa (AO) potential [eq. ^] relevant to experimental systems, (b) The competing 
interactions used in this work, /Sum + ^uy for various values of the Morse potential 
well depth /Sem as shown in the legend. The case of /Sem = is the pure Yukawa 
interaction [eq. (j3|)] contact potential of /3eY=kT and inverse Debye screening length 
Ka — 0.5. 



where Z is the colloid charge and Is is the Bjerrum length. To model the experimental 
system, we therefore fix the Debye length to an experimentally relevant value na = 0.5, 
and consider different values of the contact potential (Sey- In the experimental system 
we seek to model, van der Waals attractions are largely absent, due to solvent-colloid 
refractive index matching, the attractions are driven by the addition of non-absorbing 
polymer. 

We investigate the structure of these simulated colloidal clusters, with reference 
to the Morse clusters [35], using a new method we have developed, which we term 
the topological cluster classification (TCC) [IHIIIT]. Although we consider electrostatic 
repulsions, we expect little effect on the ground state for small m < 10 clusters [36] . 

3. Simulation and analysis 

We use a standard Brownian dynamics simulation scheme |1H]. The scheme generates 
a discrete coordinate trajectory as follows 

D ^ 

r,{t + 6t) = r,{t) + — J2 F,j{t)6t + 6rf, (5) 

where 5t is the simulation time step, D is the diffusion constant and Fjj is the pairwise 
interaction. The colloids respond to the direct interactions Fjj and the solvent-induced 
thermal fluctuations Srf are treated as a Gaussian noise with the variance given by the 
fluctuation-dissipation theorem. 

In the simulations m particles are initialised randomly at volume fraction = 
m{7ra^/6)/l^ = 0.0029 in a cubic box of side /. Periodic boundary conditions for the box 
are implemented. We study the evolution of one cluster per simulation, and consider 
the case where all particle are part of the cluster. The inter-particle interaction is 
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the Morse potential [eq. ([2])] with range parameter po = 33.06cr~^ which is truncated 
and shifted for r > 3cr, where the Morse potential is typically less than 10~^^. The 
electrostatic interactions are treated by adding a Yukawa repulsion term [eq. ([3])] for 
different values of /3eY as specified and this is also truncated and shifted for r > 3cr. This 
rather short value enables the particles to form clusters more easily starting from the 
initial randomised state. While the Yukawa repulsion has not fully decayed at r = 3a, 
the largest separation of two particle centres we are interested in is set by the size of 
the clusters. The largest separations for which we consider the effects of the Yukawa 
repulsion are the m = 13 clusters, where the maximum separation is less than 3a. Since 
the Morse interaction has no hard core, we remove the hard core component of eq. 
(|3]). Each state point is sampled with between four and twelve statistically independent 
simulation runs. 

The time-step is 6t = 0.03 simulation time units and all runs are equilibrated for 
10^ steps and run for further 10^ steps. The rather long simulation runs were required 
to be sure that, in the case of Yukawa repulsions, that the particles had ample time to 
interact with one another, and cluster. We define the Brownian time as the time taken 
for a colloid to diffuse its own radius: 



In the simulations, tb ~ 2474 time units, while in the experiments ~ 10 s [37]. The 
simulation runs therefore correspond to a total of around 68 hours, a timescale certainly 
comparable to experimental work. 

We identify two particles as bonded if the separation of the particle centres is 
less than 1.25cr, which is close to the attractive range of a strict AO potential [eq. ([T])]. 
Having identified the bond network, we use the Topological Cluster Classification (TCC) 
to determine the nature of the cluster jlTj. This analysis identifies all the shortest path 
three, four and five membered rings in the bond network. We use the TCC to find 
clusters which are global energy minima of the Morse potential. We follow Doye et.al. 
[35] and term these clusters 3A, 4A, 5A, 6A, 7A, lOB, and 13B for m = 3, 4, 5, 6, 7, 10, 13. 
For m ^ 7 there is one global minimum for the Morse potential. At higher m there 
are multiple minima corresponding to different values of the range parameter po- lOB 
and 13B correspond to a short ranged Morse potential. We assume that these are the 
relevant global minima for po = 33.06. In addition, in the case of m = 13 clusters we 
identify the FCC and HCP thirteen particle structures in terms of a central particle and 
its twelve nearest neighbours. For more details see [17] . 

4. Results 

We shall consider each size in turn, before drawing together our results. We use two 
control parameters, the well depth of the attractive Morse interaction l3eM [eq. ([2])] 
and the contact potential of the Yukawa repulsion /3ey [eq. ([3])]. Increasing f3eM thus 
promotes clustering, while /3ey suppresses clustering. In experiments on colloids, the 
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Figure 2. (color online) (a) Cluster populations for the m = 3 system as a function 
of the well depth of the attractive interaction Pem with /Sey = 0. Here we define 
cluster population by the ratio Nc/Nm where Nc and Nm are the total number of 
specific clusters and the total number of m membered clusters identified respectively 
in the simulation after equilibration. We see a clear trend towards the 3A triangle 
as the minimum energy ground state structure. The orange line corresponds to 
tl = 400tb, our criterion for ergodicity breaking on the simulation timescale [see 
(c)]. (b) population of clusters in the minimum energy state (3A triangle) for different 
strengths of the Yukawa interaction /Sey . The principal effect of the Yukawa repulsion 
is to shift the curves to require higher values of /Sem to achieve the same population 
of clusters in the ground state, (c) the average bond lifetime as a function of /Sem 
for various f^Ey as indicated. The orange lines correspond to tl = AOOtb- (d) the 
linear— >3A transition can be accomplished by rotating one end particle around the 
central particle. This process involves no energetic penalty for f^Ey — 0. 



electrostatic charge is usually not systematically varied. Therefore we consider specific 
values of f3eY and plot the response of the system to (3eM- Small clusters with m ^ 5 are 
able to reach the minimum energy ground state structures (which maximise the number 
of bonds), while for larger clusters, geometric frustration leads to kinetic trapping which 
severely limits access to the ground state at high values of the attractive interaction. 

4-1- Small clusters m ^ 5 

We begin our presentation of results by considering the m = 3 system, as shown in 
figure |2](a). The main conclusion from these data, as with all m < 5, is that increasing 
the attractive interaction strength f3eM leads to a higher population in the minimum 
energy ground state, here the '3A' triangle with D^h point group symmetry and three 



Geometric frustration in small colloidal clusters 



8 



bonds. This occurs at the expense of higher energy clusters, which, for m = 3 are 
'hnear' clusters with two bonds. In the case of Pey = 0, the potential energy for a 
short-range potential such as eq. ([2]) is approximately equal to the number of bonds. 
Since 3A triangles have three bonds, and the linear clusters have only two, there is a 
strong thermodynamic driving force to the 3A cluster as the attractions are increased, 
consistent with the behaviour seen in figure [2t^a). 

The effect of increasing the Yukawa repulsion is to slightly suppress the development 
of the 3A population [figure 121(b)], as expected in this system with competing 
interactions. In other words, the slight upwards shift in the potential (3um + Puy 
[figure [11(b)] due to the Yukawa contribution acts as a relatively small perturbation to 
the m = 3 system. However, Psy = 5 substantially suppressed the colloidal aggregation 
and few three-membered clusters were formed. We hereafter consider [iey = 0, 1, 3 only. 

The average bond lifetime tl is also shown in figure Ws). Here we define bond 
lifetime as the time between a bond formation event (where the separation between two 
colloids falls below the bond length 1.25cr) a bond breakage event (where the separation 
between the same two colloids rises above the bond length). The bond lifetimes were 
widely distributed in all cases. We see that for Psm ^ 10, is very much less than 
the simulation time, so the system may be regarded as equilibrated. In this regime, 
Tl exhibits an Arrhenius-like behaviour, as expected for an equilibrated system. At 
higher values of the interaction strength, the average bond lifetime approaches, and 
then exceeds the simulation run time, so the system is unable to equilibrate, on the 
simulation (and experimental [2Z]) timescale. Thus the system is regarded as nonergodic 
on these timescales, in that it cannot explore all its configurations. We take a average 
bond lifetime r/ = AQQtb as a crossover between ergodic and nonergodic. This ergodic 
- nonergodic transition is indicated in figure |2l(a). However, although the system is 
nonergodic, the absence of any geometric frustration enables the minimum energy 
ground state to be reached. The absence of geometric frustration is understood as 
follows [figure |2l(d)]: for m = 3, if fisy = there is no energy barrier in the linear— ;'3A 
transition, because there is no angular dependence in the interaction. In other words, a 
steepest descent quench for a three-membered cluster yields the 3A triangle. 

Turning to the case for m = 4, the situation is somewhat more complex [figure 131(a)] 
Rather than two states, there are four: 4A tetrahedra (the ground state with 6 bonds 
and T/i point group symmetry) diamonds (5 bonds), triangle-lines and squares (4 bonds) 
and linear (3 bonds). Squares are distinct from diamonds in that there are no diagonal 
bonds. However, like the m = 3 case, increasing l3eM we find a peak in the population 
of linear clusters, triangle-lines and diamonds respectively. Each higher energy state has 
a distribution which becomes progressively less favoured at higher values of f3eM, as the 
4 A becomes the dominant structure, the yield of 4A approaches unity for /3eM > 10. 
Squares have a rather low yield, much less than triangle-lines (3A-(-l), which have the 
same number of bonds. We return to the case of competing structures during our 
analysis of m = 6. 

The effect of increasing the Yukawa interaction is similar to the m = 3 case: 
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Figure 3. (color online) (a) Cluster populations for the m — A system as a function 
of the well depth of the attractive interaction (Sem with /3ey = 0. Here we see a clear 
trend towards the 4A tetrahedron as the ground state structure. We identify three 
higher energy structures: diamonds (5 bonds), triangle-lines (or 3A+1) (4 bonds) and 
linear (3 bonds). The orange line corresponds to t; = 400tb. (b) population of clusters 
in the minimum energy ground state (4A tetrahedron) and extended diamond structure 
for different strengths of the Yukawa repulsion /Sey ■ The principle effect of the Yukawa 
repulsion is to shift the curves to require higher interaction strengths to achieve the 
same degree of clusters in the ground state: the more elongated diamond structure 
does not appear more favoured for stronger Yukawa repulsions. 
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Figure 4. (color online) Cluster populations for the m = 5 system as a function of the 
well depth of the attractive interaction /Sem with /3ey = 0. The orange line corresponds 
to Ti = 400tb . sp4b denotes a four-membered ring with one particle bonded to all four 
in the ring. 



the development of the 4A population is somewhat suppressed [figure El^b)]. It has 
been suggested that the introduction of repulsions might be expected to promote more 
elongated structures [12]. For a given strength of the attractive interactions f3eM, there 
is indeed a tendency towards a higher population of the elongated 4D diamond structure, 
however the overall trend is unaltered and we restrict our analysis to the /3ey = case. 

We now consider the m = 5 system for Pey = [figure IH^a)]. We find that for 
f^^M ^ 8.0, the overwhelming majority of clusters are in the ground state, 5A (triangular 
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Figure 5. (color online) (a) Cluster populations for the m = 6 system as a function 
of the well depth of the attractive interaction /Sem with fSey = 0. Here at strong 
interaction strength, rather than the Morse ground state 6A octahedron, instead 6Z are 
found, the ground state for the Dzugutov potential, (b). The orange line corresponds 
to tl = 400tb. (c) The 3A— >-4A— >5A— )>6Z aggregation pathway with does not involve 
the breaking of any bonds and thus promotes the formation of 6Z over 6A for high 
values of ^Em where the average bond lifetime r/ 

is comparable to or greater than the simulation runtime. 

bipyramid with 9 bonds and D^h point group symmetry), in a similar way to the m = 3 
and m = 4 cases. Defective triangular bipyramids (4A+1) form an excited state with 
10 or 11 bonds. Clusters based around 4-membered rings with 5-8 bonds (sp4b and 
squares+1) are present in yields up to a few percent for relatively weak attractions 
{(3eM ~ 5). Five-membered rings (pentagons, five bonds) are found in small quantities, 
similar to the square in the case of m = 4. The main result from considering these 
small clusters is that, although the system may become non-ergodic on the simulation 
timescale, the clusters can nevertheless access their global energy minima for sufficient 
interaction strengths [Pem ^ 8.0). In other words, access to the global minimum is not 
geometrically frustrated. 

For m = 5 we do not explicitly consider all possible structures. For example linear 
clusters may form at weaker interaction strengths. At Pem = 4.5, for example, only 
around 10% of m = 5 clusters are identified, however for most interaction strengths 
we consider, the vast majority of clusters are one of the five structures considered. A 
similar argument holds for m = 6 and m = 7. 

4-2. m = 6 clusters: competing structures 

For 3 < m < 5 the dominant structure at high interaction strength is the global energy 
minimum for the Morse interaction. However, in the case of m = 6, we find this is 
not the case [figure [5]^ a)]. We see only a small population of the Morse global energy 
minimum, the 6A octahedron {Oh point group symmetry) with another structure of C2v 
point group symmetry appearing to dominate the system. This cluster is the ground 
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state for the Dzugutov potential [plotted in figure [5l^b)] [19], and we term it the 6Z. 
Why should this Dzugutov cluster be more popular than 6A octahedra? Both clusters 
have 12 bonds (near-neighours) , and for a Morse potential with a relatively short range 
as we use here, the energy contribution from second-nearest neighbours is a factor of 
5.625 X 10~^ of the total energy for the 6A (it is essentially zero for the 6Z). In fact, 
when the Yukawa repulsions are added, 6Z becomes energetically favoured as it is more 
elongated. Let us however consider the case for neutral particles where Pey = 0. 

The difference between the two clusters' ground state potential energies is slight, 
but there are two remaining contributions to the system's free energy to be considered. 
Firstly the contribution due to the vibrational modes of the clusters must be accounted 
for, and secondly, so must the volume of phase space which is accessed upon translating 
and rotating the clusters through the system volume. Assuming that the vibrational 
modes may be approximated as harmonic (valid at low temperatures), the vibrational 
contribution to the free energies can be computed using a standard normal mode 
analysis. We have done this and found the contribution due to the vibrational modes 
to be completely negligible. The contribution due to translation will be the same for 
both clusters, so need not be considered. This leaves the contribution due to rotating 
the clusters. For the rotation we must consider both the point group symmetry and the 
cluster's radius of gyration. Comparing the point group symmetries, we can see that 
there are 24 different ways to reorientate the 6A cluster which are merely a permutation 
of indistinguishable particles, while for the 6Z there are only two ways. In computing 
the entropy we must count each permutation of indistinguishable particles only once. 
For this reason, in computing the entropy, we are able to rotate the 6Z cluster through 
a greater portion of the available phase space. This results in an increase in the 6Z 
cluster's population relative to the 6A cluster's by a factor of 12. One could think of 
the 6Z cluster's lower symmetry as a form of increased disorder. We also consider the 
radius of gyration: Rq = ^^=i(i"i — i"cm)^ where are the particle coordinates and 
Tcm is the centre of mass. Rq is larger by a factor of approximately 1.06 times for the 
6Z cluster than it is for the 6A. So upon rotating the 6Z cluster its particles traverse 
a greater portion of the available phase space than is the case for the 6A cluster. This 
increases the entropy of the 6Z cluster, relative to that of the 6A, further favouring it. 
It seems plausible that this could account for the relative differences in the populations, 
i.e., that the 6Z is thermodynamically favoured over the 6A by a factor of 30 as shown 
in figure M^a). 

As PeM becomes very large, we would ultimately expect a trend towards a 6A 
dominated population due to the (small) difference in potential energy. However, on 
these simulation timescales, the average bond lifetime is too long to enable the transition 
to a 6A dominated population, and, since 6Z can be the result of a 3A^4A^5A— ;'6Z 
aggregation sequence [Fig[5]^c)], and the formation of 6A requires bond breakage, at 
strong interaction strengths, 6Z dominates for kinetic reasons. 

Like the smaller clusters, for m = 6, we identify different structures which become 
significant at lower values of (3eM- In decreasing stability, these are clusters with two 



Geometric frustration in small colloidal clusters 



12 




Figure 6. (color online) Cluster populations for the m = 7 system as a function of 
the well depth of the attractive interaction Pem with /?£y = 0. Here the Morse ground 
state 7 A pentagonal bipyramid is readily formed at intermediate interaction strength. 
The orange line corresponds to ti = AOOtb- sp5b denotes a five-membered ring with 
one particle bonded to all five in the ring. A cluster based on three overlapping 5A 
triangular bipyramids, 3x5A, shows a re-entrant behaviour, dominating the cluster 
population at high and low values of ^em- 

tetrahedra [which we denote as 5A-I-1 in figure |5]^a), with either 10 or 11 bonds], defective 
octahedra (denoted as sp4b-|-l, with between 9 and 11 bonds), defective pentagonal 
bipyramids (denoted as sp5b with 10 bonds) and clusters formed of a five-membered 
ring with one bound particle (pentagon+l, between 6 and 9 bonds). 

4.3. Larger clusters: geometric frustration 

For the small clusters we have considered so far, the number of particles is apparently 
too small to form metastable states. However, for m = 7 we see that the yield of the 
Morse global energy minimum 7A pentagonal bipyramid approaches unity for moderate 
strengths of f3eM, but for (3eM > 12 not all simulations reach the 7A [figure EJ^a)]. 
Once the 7A is formed, the system remains in a 7A state, but other metastable states 
have lifetimes longer than the simulation runs. The rise in 7A population as a function 
of Pem appears to slow around the ergodic-nonergodic transition (which we define as 

TL=mTB). 

In the nonergodic regime where bond breaking governs those structures which form, 
we expect an aggregation sequence similar to the 3A— ;'4A— )-5A— )-6Z sequence shown in 
figure M^c). Stepwise aggregation of one particle onto a 6Z cluster would lead to a 
structure we term 3x5A, as it may be decomposed into 3 overlapping 5A triangular 
bipyramids, or, equivalently, a 6Z with an additional tetrahedron. This structure has 
15 bonds. The 7A has 15 bonds where the separation ^ cr, close to the minimum of 
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Figure 7. (color online) (a) Cluster populations for the m = 10 (a) and w = 13 (b) 
system as a function of the well depth of the attractive interaction Pem with /^ey — 0. 
The orange line in (a) corresponds to t/ — 400ts 



the Morse potential, while the distance between the top and bottom particle is around 
1.02283(7, contributing 0.7193£m to the energy. The 7A is thus around 0.7193/3eAf lower 
in energy than the 3x5A structure. In the nonergodic regime for ^em > 12, we see 
that this 3 x 5A structure dominates the system, due to kinetic trapping. The behaviour 
of the 7A system in the nonergodic regime highlights the sampling limitations of our 
simulation approach. To map this regime more accurately, many more simulations are 
desirable than the 4 runs per state point we have been able to perform. 

At weaker interaction strengths, like the smaller clusters, a variety of structures 
become popular, based on a diminishing number of bonds. Of these, sp5b+l, a 
defective pentagonal bipyramid with 11-13 bonds shows a rather slow rate of decay 
upon increasing /3eM- At the weakest interaction strength, f3eM = 5.5, the yield of 3x5A 
exceeds that of 7A. This 're-entrant' behaviour of 3x5A is apparently a consequence of 
the different number of bonds relative to 7A. Cluster with fewer bonds are promoted 
for weak interaction strengths, but in the nonergodic regime, kinetic trapping promotes 
3x5A. 

At larger sizes again we consider the Morse global minima, which are lOB for 13B 
for m = 10, 13 clusters respectively. Once more the effects of frustration are clear. The 
m = 10 system (figure [7j) features a clear maximum yield of lOB at Pem ~ 7, however 
the yield is only around 14%. The rise in the yield of lOB as a function of Pem up to 
the ergodic-nonergodic crossover suggests that this rise in average bond lifetime is the 
primary mechanism by which a further rise in the lOB yield is suppressed. In other 
words, kinetic trapping prevents access to lOB at higher strengths of the attractive 
interaction. 

For m = 13, in addition to 13B, we also find crystal fragments. The maximum yield 
of 13B is around 10%. For I3em = 7.0, HOP crystal fragments are in fact more popular 
than the 13B ground state. We find no icosahedra, these are the global energy minimum 
for longer ranged Morse interactions (po < 14.76) [35] than we use here, although for 
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Figure 8. (color online) (a) Elongation d/Rc for the m = 10 (a) and m = 13 (b) 
system as a function of the well depth of the attractive interaction Pem with /^ey = 0. 
The dashed orange line denotes the elongation in the case of the minimum energy 
ground state . 



moderate values of ^Sm it is not unreasonable to expect icosahedra. We expect that a 
similar kinetic argument to that suggested above for m = 10 concerning the low yield of 
lOB holds for m = 13 as well. For m = 10 and m = 13 we cannot exclude the possibility 
that other global minima exist. Doye et. al. [35] calculated global minima for po ^ 25. 
The smaller global energy minima for m ^ 6 exhibit no strain for short-ranged Morse 
interactions, and only a limited amount of strain in the case of 7A, these structures are 
therefore also expected to be minima for po = 33.06 as used here. This not necessarily 
the case for m = 10 and m = 13. We are however unaware of any more appropriate 
structures than those listed in [35] . 

A question arises in comparing figure [7] with results for smaller clusters, such as 
the case of m = 3 (figure [2]). In general one expects that larger clusters should be able 
to form at higher temperatures (weaker interaction strengths) [1], (although for small 
Lennard- Jones clusters, the melting line is in fact non-monotonic as a function of m 
[50]). Here, we are interested in whether all m particles in the simulation box aggregate 
to form a cluster. This is less likely for larger m, so our statistics suffer for lower values 
of relative to the smaller clusters considered. It is nevertheless very clear that the 
yield of the assumed global minimum clusters is markedly reduced in the case of m = 10 
and m = 13 relative to smaller clusters. We now proceed to consider elongation. 



4.4- Elongation 

The ground states of larger clusters have been found to be strongly affected by the 
strength of a Yukawa repulsion, with stronger repulsions leading to more elongated 
structures [36]. This effective long-ranged repulsion was indeed predicted to have a 
profound effect upon both the shape and the size of clusters of charged colloids by 
Groenewold and Kegel [12] , and experiments have found evidence for elongated clusters 
[SI E] and Bernal spirals [0]. We therefore investigate the degree of elongation of the 
larger clusters considered here. We consider the largest separation of two coordinates 
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within a cluster, d, and divide this by the radius of gyration Rg of the same coordinate 
set. Thus, larger d/Rc corresponds to more elongated clusters, and for a large spherical 
cluster with m — )■ oo, d/Rc 2 a/5/3. For highly symmetric clusters such as icosahedra, 
d/Rc = 2.08. 

Figure M shows the elongation d/Rc as a function of f3eM- The elongation for 
m = 10 and m = 13 appears to be dominated by failure to access the ground state; 
nonergodic systems are more elongated. The data for /3ey = 1.0 are included in figure 
IHl while for /3ey = 3.0, few clusters formed at all, and none in the ergodic regime. From 
the discussion above, one might expect the /3£y = 1 data to show a higher degree of 
elongation due to the Yukawa repulsion, however figure M^a) indicates that this is not 
the case for m = 10 and m = 13, except for a slight hint at the lower values of Pem- 
We might expect that elongation induced by long-ranged repulsions is promoted by 
increasing m, or Pey or perhaps by a longer-ranged attraction that might suppress the 
transition to non-ergodicity, due to a less complex potential energy surface (see section 
E]) [51]. However in this system we see little evidence for elongation. 

5. Discussion 

The behaviour of the systems considered here can readily be decomposed into systems 
where frustration is not relevant, m ^ 5 where a well-defined structure is favoured, and 
those where frustration leads to a non-trivial potential energy landscape and limited 
access to the ground state m ^ 7. In the case of m = 6 two structures with the same 
number of bonds, 6A octahedra and 6Z D^h compete. Of particular relevance here is 
the short range of the attractive interactions. It has been noted previously that these 
tend to promote a complex energy landscape, for example m = 13 Morse clusters with 
po = 14.0 have some 54439 local energy minima, compared to just 685 for the longer- 
ranged case of Po = 4.0 [51j. The Morse potential with po = 14.0 is an approximation to 
Ceo, clusters of which are known to exhibit kinetic trapping [52], although in that case, 
icosahedral clusters were favoured. With po = 33.06, stronger trapping is likely. What 
is clear, therefore, is that these systems can, only in the simplest cases (m = 3,4,5), 
form high yields of clusters in the minimum energy ground state, although m = 7 does 
have a window in (3em where the yield of 7A is rather substantial. 

5.1. Relevance to experiments 

This work has been largely motivated by experimental studies of clustering in colloidal 
dispersions. In those systems, prized for their tunable interactions, the attractions, 
either Asakura-Oosawa, or van der Waals are typically rather short ranged, and thus 
likely to exhibit kinetic trapping similar to the systems studied here. We therefore argue 
that for substantial yields of more complex 'colloidal molecules', it is appropriate to seek 
more sophisticated means than the spherically symmetric spheres we have considered 
here, such as patchy particles [101 HH [13 SIl ES] and Janus particles [ISl [IZ], or to 
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design routes of preparation |ll|5lj. However, even 'purpose-designed' patchy particles 
can have only a rather limited window where the yield of ground state structures is 
substantial [Tn[T2|[T3]. 

We now compare our results to those found in an experimental system in which 
the interactions are, to a first approximation, identical Each of the simulations 

was conducted with a fixed number of of particles corresponding to the cluster size 
investigated. By contrast, most experiments are carried out with a bulk colloidal 
suspension, in which clusters of different size form. The experimental data report the 
relative abundance of cluster types for a fixed number of particles per cluster, m, and 
each cluster is assumed to interact only weakly with other clusters. In the experimental 
system, m = 3 formed a majority of clusters in the 3A triangle which maximises the 
number of bonds. This occurred at interaction strengths comparable to those found in 
the simulations reported here. In the experiments m = 4 and 5 formed only around 
10% of 4A tetrahedra and 5A triangular bipyramids respectively, in sharp contrast to 
the simulation results, where the yield was essentially 100%. In the case of m = 6, 
both in experiment and simulation, 6Z C2v dominated the 6A octahedron. However, in 
the simulation, in the ergodic regime, the population difference was around a factor of 
30, while in the experiment the difference was at least an order of magnitude larger. 
Furthermore, the maximum yield of 6Z was an order of magnitude lower in the 
experiment than the simulation, where it was around 100%. For m = 7 both experiment 
and simulation show signs of geometric frustration, with the yield of 7A being reduced 
in the nonergodic regime of strong interaction strength. However, while the peak yield 
of 7A approached 100% in the simulation, the experiments were limited to a few percent 
of 7A. We consider possible origins for this discrepancy in the following section. 

5.2. Charging in apolar colloidal systems 

Before concluding, it is worth considering these results in the light of some other recent 
experimental studies. Campbell et. al. [B] report clusters in a system in which they 
measured the colloid charge in a dilute suspension to be Z = 140 e per 1.5 /im diameter 
colloid, where e is the elementary charge. According to eq. (ED, this maps to a Yukawa 
contact potential Pey = 35, their quoted value for the Debye length is comparable to 
ours. We have observed only very limited clustering at Psy = 5, corresponding to a 
charge of Z = 47 e and expect none at higher Yukawa contact potentials. Dibble et. 
al. [8J quote a similar value of Z = 165 e per 2 /im colloid. Moreover Sedgwick et. al. [7] 
report Z < 10^ e in their study of clustering. Although not strictly inconsistent, this 
seems rather higher than the values we predict from eq. i^. 

In their simulation study of gelation in colloidal systems with competing 
interactions, Sciortino et. al. |55] used a comparable contact potential for the Yukawa 
repulsion Pey to that used in the simulations here. In other words, a much weaker 
Yukawa repulsion than that expected from the colloid charge quoted by Campbell 
et. al. [6]. Interestingly, Sciortino et. al. f55] found that similar Bernal spiral 
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structures were formed in their simulations to those observed experimentally p]. From 
the arguments presented above, it might be supposed that the Bernal spiral can 
form without significant geometric frustration. An important question remains in 
the role of colloid concentration. We have considered a rather low volume fraction 
= m{7Ta^/6)/l^ = 0.0029, as we are interested in isolated clusters. Whether the 
details of our results apply at finite concentration where clusters interact with one 
another remains to be seen, but many colloid volume fractions quoted in the literature 
are < 0.2. In this regime we expect the discussion above regarding clustering in the 
presence of electrostatic repulsions to be reasonable. 

Given the success of linear Poisson-Boltzmann theory in describing electrostatic 
interactions in these systems in the absence of polymer-induced attractions [56], one 
is tempted to enquire as to the origin of the discrepancy between the findings of this 
paper and the experimental literature P, [TJ IH]. As far as we know, in only one case has 
quantitative agreement been found between experiment and simulation of competing 
interactions using conventional models of electrostatic repulsion and polymer-induced 
attraction [57], and there the electrostatic repulsions were screened by salt. Combined 
with the discrepancies presented here and those in the experimental system [37| to which 
we have tuned the interactions, we believe that there is more than meets the eye to these 
colloidal model systems with competing interactions. We see little reason to suppose 
that the polymer-induced attractions deviate substantially from theory therefore 
we speculate that the electrostatic repulsions in the clusters may deviate from those 
deduced by electrophoretic mobility measurements in dilute suspensions [HI El EZl ES] ■ 

6. Conclusions 

We have studied isolated colloidal clusters using Brownian dynamics simulations. Our 
system is tuned to match experimental work in colloidal systems with polymer-induced 
depletion attractions. For sufficient strengths of the attractive interaction, the average 
bond lifetime exceeds the simulation runtime, and the system as nonergodic on this 
timescale; this conclusion also applies to some experimental systems. However, for 
small clusters m ^ 5 this ergodicity breaking does not prevent the system reaching 
the minimum energy ground state structure, as no bonds need be broken. Thus for 
these small clusters, we can direct the system in a controlled way towards a prescribed 
ground state. For m = 6 we find a structure which appears to minimise the free energy 
with symmetry, rather than the octahedron that forms the minimum energy ground 
state [351 ES]- The energy difference between these structures is negligible, and for 
moderate interaction strengths (the ergodic regime) the structure is favoured for 
entropic reasons, while in the nonergodic regime, it is kinetically favoured as it is the 
product of an aggregation sequence that does not involve bond breakage. At larger 
cluster sizes, the system is kinetically frustrated from reaching the ground state in all 
but a limited window. Moreover, recent experimental results in a similar system show 
much lower yields of clusters in the ground state than we find here [37] . The origin of the 
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discrepancy is likely related to a difference in the effective interaction potentials between 
experiment and simulation. It is perhaps surprising that such small systems exhibit this 
degree of kinetic trapping; only those sizes for which there is no geometric frustration 
(m = 3,4,5) reach the minimum energy ground state. These results suggest that high 
yields of complex microstructures of colloidal clusters and molecules may benefit from 
reversible quenching [38[ \39\ SOj HI], rather than the fixed interactions of many colloidal 
systems, which lead to 'instantaneous quenching'. 

Some pointers for further work are considered. We have employed a simple approach 
in our simulations, as we are motivated to reproduce the recent experimental system. 
One approach to investigate larger clusters could be to run a simulation of a much 
larger number of particles, and to consider each cluster as a separate system. In this 
way, larger clusters could be accessed than was the case here. Although this could 
in principle resolve some discrepancies between the results presented here and those 
reported for the experimental system, we believe this is a most unlikely scenario. Other 
possibilities might be to implement the methods found in the (atomic) cluster literature 
[11 |35l |50] to comprehensively determine the structure of larger clusters. Another 
possibility would be to determine the phase diagram for the system considered here, 
using normal mode analysis to provide a theoretical prediction of the population levels 
for the various structures considered. 

To provide further support to experimentalists in the quest to control the structure 
of colloidal clusters it might be helpful to move beyond the one-component description 
employed here. Given that the charging number quoted in some experimental work is 
so small {Z < 100 e)[6l [3 El EZl [58], it might be possible to determine structures of 
colloidal clusters by for example developing the primitive model of electrolyte systems 
to colloidal systems with many ionisable sites on each particle [59]. 
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